module step_mod

  use const_mod
	use namelist_mod
	use mesh_mod
	use state_mod
	use tracer_mod
	use parallel_mod
	use time_mod, only: old=>old_time_idx, new=>new_time_idx
	use vert_coord_mod
	use operators_mod
	use math_mod
	use upwind_mod
	use ffsl_mod

	implicit none

	private

	public step_init
	public step_final
	public step_adv_split
	public step_adv_nosplit
	public step_adv_ieva

	integer :: hori_rk_nstep, vert_rk_nstep 
	logical, allocatable :: do_ieva_vert(:)

	interface
		subroutine calc_tracer_mass_flux_hori_interface(tracer, state, old_q, rk_q, qmf_lon, qmf_lat, rk_substep, do_limiter)
			import tracer_type, state_type, r8, mesh
			type(tracer_type), intent(inout) :: tracer
		  type(state_type), intent(in) :: state
		  real(r8), intent(in) :: old_q(mesh%full_lon_lb:mesh%full_lon_ub, &
		  												      mesh%full_lat_lb:mesh%full_lat_ub, &
		  												      mesh%full_lev_lb:mesh%full_lev_ub)
		  real(r8), intent(in) :: rk_q(mesh%full_lon_lb:mesh%full_lon_ub, &
		  												     mesh%full_lat_lb:mesh%full_lat_ub, &
		  												     mesh%full_lev_lb:mesh%full_lev_ub)
		  real(r8), intent(out) :: qmf_lon(mesh%half_lon_lb:mesh%half_lon_ub, &
		  														     mesh%full_lat_lb:mesh%full_lat_ub, &
		  														     mesh%full_lev_lb:mesh%full_lev_ub)
		  real(r8), intent(out) :: qmf_lat(mesh%full_lon_lb:mesh%full_lon_ub, &
		  														     mesh%half_lat_lb:mesh%half_lat_ub, &
		  														     mesh%full_lev_lb:mesh%full_lev_ub)
		  real(8), intent(in) :: rk_substep
		  logical, intent(in), optional :: do_limiter
		end subroutine calc_tracer_mass_flux_hori_interface

		subroutine calc_tracer_mass_flux_vert_interface(tracer, state, old_q, rk_q, we, qmf_lev, rk_substep, do_limiter)
			import tracer_type, state_type, r8, mesh
			type(tracer_type), intent(inout) :: tracer
		  type(state_type), intent(inout) :: state
			real(r8), intent(in) :: old_q(mesh%full_lon_lb:mesh%full_lon_ub, &
			 	 			  						        mesh%full_lat_lb:mesh%full_lat_ub, &
			 	 			  				  		      mesh%full_lev_lb:mesh%full_lev_ub)
		  real(r8), intent(inout) :: rk_q(mesh%full_lon_lb:mesh%full_lon_ub, &
			 	 			  								      mesh%full_lat_lb:mesh%full_lat_ub, &
			 	 			  								      mesh%full_lev_lb:mesh%full_lev_ub)
		  real(r8), intent(in) :: we(mesh%full_lon_lb:mesh%full_lon_ub, &
															   mesh%full_lat_lb:mesh%full_lat_ub, &
															   mesh%half_lev_lb:mesh%half_lev_ub)
		  real(r8), intent(out) :: qmf_lev(mesh%full_lon_lb:mesh%full_lon_ub, &
																       mesh%full_lat_lb:mesh%full_lat_ub, &
																       mesh%half_lev_lb:mesh%half_lev_ub)
		  real(8), intent(in) :: rk_substep
		  logical, intent(in), optional :: do_limiter

		end subroutine calc_tracer_mass_flux_vert_interface
	end interface

	procedure(calc_tracer_mass_flux_hori_interface), pointer :: calc_tracer_mass_flux_hori => NULL()
 	procedure(calc_tracer_mass_flux_vert_interface), pointer :: calc_tracer_mass_flux_vert => NULL()

contains
	
	subroutine step_init()

		if (adv_hori_scheme == 'upwind') then
			hori_rk_nstep = 3
			calc_tracer_mass_flux_hori => calc_tracer_mass_flux_hori_upwind
		else if (adv_hori_scheme == 'ffsl') then
			hori_rk_nstep = 1
			calc_tracer_mass_flux_hori => calc_tracer_mass_flux_hori_ffsl
		end if

		if (adv_vert_scheme == 'upwind') then
			vert_rk_nstep = 3
			calc_tracer_mass_flux_vert => calc_tracer_mass_flux_vert_upwind
	  else if (adv_vert_scheme == 'ffsl') then
			vert_rk_nstep = 1
			calc_tracer_mass_flux_vert => calc_tracer_mass_flux_vert_ffsl
		end if 

		if (.not. allocated(do_ieva_vert)) then
			allocate(do_ieva_vert(3))
			do_ieva_vert(1:3) = .false.
			do_ieva_vert(3) = tracer_vert_ieva
		end if

		call ffsl_init()
		call upwind_init()

	end subroutine step_init

	subroutine step_final()

		if (allocated(do_ieva_vert)) deallocate(do_ieva_vert)

		call upwind_final()

	end subroutine step_final

	subroutine step_adv_nosplit(old_state, new_state, old_tracer, new_tracer, dt)

		type(state_type), intent(inout) :: old_state
		type(state_type), intent(inout) :: new_state
		type(tracer_type), intent(inout) :: old_tracer
		type(tracer_type), intent(inout) :: new_tracer
		real(r8), intent(in) :: dt
		real(r8) rk_q(mesh%full_lon_lb:mesh%full_lon_ub, &
								  mesh%full_lat_lb:mesh%full_lat_ub, &
								  mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8) pole, rk_substep
		integer i, j, k, itracer, irk_step
		integer, parameter :: nrk_step= 3

		call calc_hori_uv(old_tracer, old_state, new_state, dt)
		if (mesh%num_full_lev > 1) then
			call calc_ph(old_state)
    	call calc_m(old_state)
			call calc_vert_we(old_tracer, old_state, new_state)
		else
			old_state%m = 1.0
		end if
		new_state%m = old_state%m

		associate (old_q  => old_tracer%q,       &
							 old_m  => old_state%m,        &
							 new_m  => new_state%m,        &
							 we     => old_tracer%we     , &
			        qmf_lon => old_tracer%qmf_lon, &
			        qmf_lat => old_tracer%qmf_lat, &
			        qmf_lev => old_tracer%qmf_lev)
		
		do itracer = 1, size(old_tracer%q, 4)
			rk_q = old_q(:,:,:,itracer)
			do irk_step = 1, nrk_step
	 			rk_substep = dt / (nrk_step + 1 - irk_step)
	 			if (adv_hori_scheme == 'upwind') then
			  	call calc_tracer_mass_flux_hori_upwind(old_tracer, old_state, old_q(:,:,:,itracer), rk_q, qmf_lon, qmf_lat, rk_substep, do_hori_limiter(irk_step))
			  else if (adv_hori_scheme == 'ffsl') then
			  	call calc_tracer_mass_flux_hori_ffsl(old_tracer, old_state, old_q(:,:,:,itracer), rk_q, qmf_lon, qmf_lat, rk_substep)
			  end if
		  	if (adv_vert_scheme == 'upwind') then
			  	call calc_tracer_mass_flux_vert_upwind(old_tracer, old_state, old_q(:,:,:,itracer), rk_q, we, qmf_lev, rk_substep, do_vert_limiter(irk_step))
				else if (adv_vert_scheme == 'ffsl') then
					call calc_tracer_mass_flux_vert_ffsl(old_tracer, old_state, old_q(:,:,:,itracer), rk_q, we, qmf_lev, rk_substep)
				end if

	 			do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				  do j = mesh%full_lat_ibeg + 1, mesh%full_lat_iend - 1
				  	do i = mesh%full_lon_ibeg, mesh%full_lon_iend
				  		rk_q(i,j,k) = (old_m(i,j,k) * old_q(i,j,k,itracer) - (&
				  								(qmf_lon(i,j,k) - qmf_lon(i-1,j,k)) * old_m(i,j,k) * mesh%le_lon(j)  + &
				  								(qmf_lat(i,j  ,k) * mesh%le_lat(j  ) - & 
				  								 qmf_lat(i,j-1,k) * mesh%le_lat(j-1)) * old_m(i,j,k)) / mesh%area_cell(j) * rk_substep - &
				  								 (qmf_lev(i,j,k+1) - qmf_lev(i,j,k))) / new_m(i,j,k)
				  	end do
				  end do
				end do

				! South pole
			  do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				  j = mesh%full_lat_ibeg
				  pole = sum(qmf_lat(mesh%full_lon_ibeg:mesh%full_lon_iend,j,k)) * mesh%le_lat(j) / mesh%num_full_lon / mesh%area_cell(j)
				  do i = mesh%full_lon_ibeg, mesh%full_lon_iend
				  	rk_q(i,j,k) = (old_m(i,j,k) * old_q(i,j,k,itracer) - pole * old_m(i,j,k) * rk_substep - &
				  	              								 (qmf_lev(i,j,k+1) - qmf_lev(i,j,k))) / new_m(i,j,k)
				  end do
				end do
			  ! North pole
			  do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				  j = mesh%full_lat_iend
				  pole = sum(qmf_lat(mesh%full_lon_ibeg:mesh%full_lon_iend,j-1,k)) * mesh%le_lat(j-1) / mesh%num_full_lon / mesh%area_cell(j)
				  do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						rk_q(i,j,k) = (old_m(i,j,k) * old_q(i,j,k,itracer) + pole * old_m(i,j,k) * rk_substep - &
						              								 (qmf_lev(i,j,k+1) - qmf_lev(i,j,k))) / new_m(i,j,k)
				  end do
				end do
				call fill_zonal_halo_cell(rk_q, all_halo=.true.)
				call fill_merid_halo_cell(rk_q, all_halo=.true.)
			end do

			do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				do j = mesh%full_lat_ibeg, mesh%full_lat_iend
					do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						new_tracer%q(i,j,k,itracer) = rk_q(i,j,k)
					end do
				end do
			end do
			call fill_merid_halo_cell(new_tracer%q(:,:,:,itracer), all_halo=.true.)
			call fill_zonal_halo_cell(new_tracer%q(:,:,:,itracer), all_halo=.true.)
	  end do
		end associate
	end subroutine step_adv_nosplit

	subroutine step_adv_split(old_state, new_state, old_tracer, new_tracer, dt)

		type(state_type), intent(inout) :: old_state
		type(state_type), intent(inout) :: new_state
		type(tracer_type), intent(inout) :: old_tracer
		type(tracer_type), intent(inout) :: new_tracer
		real(r8), intent(in) :: dt
		real(r8) h_q(mesh%full_lon_lb:mesh%full_lon_ub, &
								 mesh%full_lat_lb:mesh%full_lat_ub, &
								 mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8) pole, rk_substep
		integer i, j, k, itracer, irk_step

		call calc_hori_uv(old_tracer, old_state, new_state, dt)

		if (mesh%num_full_lev > 1) then
			call calc_ph(old_state) 
			call calc_m(old_state)
			call calc_vert_we(old_tracer, old_state, new_state)
		else
			old_state%m = 1.0
		end if
		new_state%m = old_state%m

		associate (old_q  => old_tracer%q,       &
							 old_m  => old_state%m,        &
							 new_m  => new_state%m,        &
							 we     => old_tracer%we,      &
			        qmf_lon => old_tracer%qmf_lon, &
			        qmf_lat => old_tracer%qmf_lat, &
			        qmf_lev => old_tracer%qmf_lev)

		do itracer = 1, size(old_tracer%q, 4)
			h_q = old_q(:,:,:,itracer)
			do irk_step = 1, hori_rk_nstep
	 			rk_substep = dt / (hori_rk_nstep + 1 - irk_step)
	 			if (adv_hori_scheme == 'upwind') then
			  	call calc_tracer_mass_flux_hori_upwind(old_tracer, old_state, old_q(:,:,:,itracer), h_q, qmf_lon, qmf_lat, rk_substep, do_hori_limiter(irk_step))
			  else if (adv_hori_scheme == 'ffsl') then
			  	call calc_tracer_mass_flux_hori_ffsl(old_tracer, old_state, old_q(:,:,:,itracer), h_q, qmf_lon, qmf_lat, rk_substep)
			  end if

	 			do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				  do j = mesh%full_lat_ibeg + 1, mesh%full_lat_iend - 1
				  	do i = mesh%full_lon_ibeg, mesh%full_lon_iend
				  		h_q(i,j,k) = (old_m(i,j,k) * old_q(i,j,k,itracer) - (&
				  								(qmf_lon(i,j,k) - qmf_lon(i-1,j,k)) * old_m(i,j,k) * mesh%le_lon(j)  + &
				  								(qmf_lat(i,j  ,k) * mesh%le_lat(j  ) - & 
				  								 qmf_lat(i,j-1,k) * mesh%le_lat(j-1)) * old_m(i,j,k)) / mesh%area_cell(j) * rk_substep) / new_m(i,j,k)
				  	end do
				  end do
				end do
				! South pole
			  do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				  j = mesh%full_lat_ibeg
				  pole = sum(qmf_lat(mesh%full_lon_ibeg:mesh%full_lon_iend,j,k)) * mesh%le_lat(j) / mesh%num_full_lon / mesh%area_cell(j)
				  do i = mesh%full_lon_ibeg, mesh%full_lon_iend
				  	h_q(i,j,k) = (old_m(i,j,k) * old_q(i,j,k,itracer) - pole * old_m(i,j,k) * rk_substep) / new_m(i,j,k)
				  end do
				end do
			  ! North pole
			  do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				  j = mesh%full_lat_iend
				  pole = sum(qmf_lat(mesh%full_lon_ibeg:mesh%full_lon_iend,j-1,k)) * mesh%le_lat(j-1) / mesh%num_full_lon / mesh%area_cell(j)
				  do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						h_q(i,j,k) = (old_m(i,j,k) * old_q(i,j,k,itracer) + pole * old_m(i,j,k) * rk_substep) / new_m(i,j,k)
				  end do
				end do
				call fill_zonal_halo_cell(h_q, all_halo=.true.)
				call fill_merid_halo_cell(h_q, all_halo=.true.)
			end do

			if (mesh%num_full_lev > 1) then
				do k = mesh%full_lev_ibeg, mesh%full_lev_iend
					do j = mesh%full_lat_ibeg, mesh%full_lat_iend
						do i = mesh%full_lon_ibeg, mesh%full_lon_iend
				  		old_q(i,j,k,itracer) = h_q(i,j,k)
				  	end do
				  end do
			  end do

	 			do irk_step = 1, vert_rk_nstep
	 				rk_substep = dt / (vert_rk_nstep + 1 - irk_step)
	 				if (adv_vert_scheme == 'upwind') then
				  	call calc_tracer_mass_flux_vert_upwind(old_tracer, old_state, old_q(:,:,:,itracer), h_q, we, qmf_lev, rk_substep, do_vert_limiter(irk_step))
					else if (adv_vert_scheme == 'ffsl') then
						call calc_tracer_mass_flux_vert_ffsl(old_tracer, old_state, old_q(:,:,:,itracer), h_q, we, qmf_lev, rk_substep)
					end if
					do k = mesh%full_lev_ibeg, mesh%full_lev_iend
						do j = mesh%full_lat_ibeg, mesh%full_lat_iend
							do i = mesh%full_lon_ibeg, mesh%full_lon_iend
								h_q(i,j,k) = (old_m(i,j,k) * old_q(i,j,k,itracer) - &
														  (qmf_lev(i,j,k+1) - qmf_lev(i,j,k))) / new_m(i,j,k)
							  if (abs(h_q(i,j,k)) > 10) then 
									print*, 'q exceeds 10'
									stop
								end if 
							end do
						end do
			    end do
			  end do
			end if 

			do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				do j = mesh%full_lat_ibeg, mesh%full_lat_iend
					do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						new_tracer%q(i,j,k,itracer) = h_q(i,j,k)
					end do
				end do
			end do
		  call fill_merid_halo_cell(new_tracer%q(:,:,:,itracer), all_halo=.true.)
			call fill_zonal_halo_cell(new_tracer%q(:,:,:,itracer), all_halo=.true.)
	  end do
		end associate
	end subroutine step_adv_split

	subroutine step_adv_ieva(old_state, new_state, old_tracer, new_tracer, dt)

		type(state_type), intent(inout) :: old_state
		type(state_type), intent(inout) :: new_state
		type(tracer_type), intent(inout) :: old_tracer
		type(tracer_type), intent(inout) :: new_tracer
		real(8), intent(in) :: dt
		real(r8) trid_a(mesh%num_full_lev), trid_b(mesh%num_full_lev), &
						 trid_c(mesh%num_full_lev), trid_r(mesh%num_full_lev), &
						 trid_out(mesh%num_full_lev)
	  real(r8) mq(mesh%full_lon_ibeg:mesh%full_lon_iend, &
	  						mesh%full_lat_ibeg:mesh%full_lat_iend, &
	  						mesh%full_lev_ibeg:mesh%full_lev_iend)
		real(r8) rk_q(mesh%full_lon_lb:mesh%full_lon_ub, &
								  mesh%full_lat_lb:mesh%full_lat_ub, &
								  mesh%full_lev_lb:mesh%full_lev_ub)	  
	  real(r8) pole, rk_substep
		integer i, j, k, itracer, irk_step

		call calc_hori_uv(old_tracer, old_state, new_state, dt)
		if (mesh%num_full_lev > 1) then
			call calc_ph(old_state)
    	call calc_m(old_state)
			call calc_vert_we(old_tracer, old_state, new_state)
		else
			old_state%m = 1.0
		end if
		new_state%m = old_state%m

		associate (old_m  => old_state%m       , &
							 new_m  => new_state%m       , &
							 old_q  => old_tracer%q      , &
							 qmf_lon=> old_tracer%qmf_lon, &
							 qmf_lat=> old_tracer%qmf_lat, &
							 qmf_lev=> old_tracer%qmf_lev, &
							 we     => old_tracer%we     , &
							 we_imp => old_tracer%we_imp , &
							 we_exp => old_tracer%we_exp )

		do itracer = 1, size(old_tracer%q, 4)
			rk_q = old_q(:,:,:,itracer)
			do irk_step = 1, hori_rk_nstep
	 			rk_substep = dt / (hori_rk_nstep + 1 - irk_step)
	 	  	call calc_tracer_mass_flux_hori(old_tracer, old_state, old_q(:,:,:,itracer), rk_q, qmf_lon, qmf_lat, rk_substep, do_hori_limiter(irk_step))

	 			do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				  do j = mesh%full_lat_ibeg + 1, mesh%full_lat_iend - 1
				  	do i = mesh%full_lon_ibeg, mesh%full_lon_iend
				  		mq(i,j,k) = old_m(i,j,k) * old_q(i,j,k,itracer) - (&
				  								(qmf_lon(i,j,k) - qmf_lon(i-1,j,k)) * old_m(i,j,k) * mesh%le_lon(j)  + &
				  								(qmf_lat(i,j  ,k) * mesh%le_lat(j  ) - & 
				  								 qmf_lat(i,j-1,k) * mesh%le_lat(j-1)) * old_m(i,j,k)) / mesh%area_cell(j) * rk_substep
				  	end do
				  end do
				end do

				!south pole
			  do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				  j = mesh%full_lat_ibeg
				  pole = sum(qmf_lat(mesh%full_lon_ibeg:mesh%full_lon_iend,j,k)) * mesh%le_lat(j) / mesh%num_full_lon / mesh%area_cell(j)
				  do i = mesh%full_lon_ibeg, mesh%full_lon_iend
				  	mq(i,j,k) = old_m(i,j,k) * old_q(i,j,k,itracer) - pole * old_m(i,j,k) * rk_substep
				  end do
				end do
			  ! North pole
			  do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				  j = mesh%full_lat_iend
				  pole = sum(qmf_lat(mesh%full_lon_ibeg:mesh%full_lon_iend,j-1,k)) * mesh%le_lat(j-1) / mesh%num_full_lon / mesh%area_cell(j)
				  do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						mq(i,j,k) = old_m(i,j,k) * old_q(i,j,k,itracer) + pole * old_m(i,j,k) * rk_substep
				  end do
				end do

				do k = mesh%full_lev_ibeg, mesh%full_lev_iend
					do j = mesh%full_lat_ibeg, mesh%full_lat_iend
						do i = mesh%full_lon_ibeg, mesh%full_lon_iend
							rk_q(i,j,k) = mq(i,j,k) / new_m(i,j,k)
						end do
					end do
				end do
				call fill_zonal_halo_cell(rk_q, all_halo=.true.)
				call fill_merid_halo_cell(rk_q, all_halo=.true.)
			end do

			if (mesh%num_full_lev > 1) then
				do k = mesh%full_lev_ibeg, mesh%full_lev_iend
					do j = mesh%full_lat_ibeg, mesh%full_lat_iend
						do i = mesh%full_lon_ibeg, mesh%full_lon_iend
							old_q(i,j,k,itracer) = rk_q(i,j,k)
						end do
					end do
				end do

	 			do irk_step = 1, vert_rk_nstep
	 				rk_substep = dt / (vert_rk_nstep + 1 - irk_step)
				  we_exp = we
				  if (do_ieva_vert(irk_step)) then
				  	call calc_we_split_wgt(old_tracer, old_state, rk_substep)
				  end if
				  call calc_tracer_mass_flux_vert(old_tracer, old_state, old_q(:,:,:,itracer), rk_q, we_exp, qmf_lev, rk_substep, do_vert_limiter(irk_step))
					do k = mesh%full_lev_ibeg, mesh%full_lev_iend
						do j = mesh%full_lat_ibeg, mesh%full_lat_iend
							do i = mesh%full_lon_ibeg, mesh%full_lon_iend
								rk_q(i,j,k) = (mq(i,j,k) - (qmf_lev(i,j,k+1) - qmf_lev(i,j,k))) / old_m(i,j,k)
							end do
						end do
					end do

					if (do_ieva_vert(irk_step)) then
						do j = mesh%full_lat_ibeg, mesh%full_lat_iend
						  do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						  	trid_a(1) = 0.0
						  	trid_c(mesh%full_lev_iend) = 0.0
							  do k = mesh%full_lev_ibeg + 1, mesh%full_lev_iend - 1
									trid_a(k)   = - rk_substep * we_imp(i,j,k) * (1.0 + sign(1._r8, we_imp(i,j,k)))
									trid_c(k-1) =   rk_substep * we_imp(i,j,k) * (1.0 - sign(1._r8, we_imp(i,j,k)))
								end do
								do k = mesh%full_lev_ibeg, mesh%full_lev_iend
									trid_b(k) = rk_substep * we_imp(i,j,k+1) * (1.0 + sign(1._r8, we_imp(i,j,k+1))) - &
															rk_substep * we_imp(i,j,k  ) * (1.0 - sign(1._r8, we_imp(i,j,k  ))) + 2.0 * old_m(i,j,k)
									trid_r(k) = 2.0 * rk_q(i,j,k) * old_m(i,j,k)
						 	  end do
								call solve_tridiagonal(trid_a, trid_b, trid_c, trid_r, trid_out, mesh%num_full_lev)
								rk_q(i,j,mesh%full_lev_ibeg:mesh%full_lev_iend)= trid_out
							end do
						end do
					end if
				end do

				do k = mesh%full_lev_ibeg, mesh%full_lev_iend
					do j = mesh%full_lat_ibeg, mesh%full_lat_iend
						do i = mesh%full_lon_ibeg, mesh%full_lon_iend
							new_tracer%q(i,j,k,itracer) = rk_q(i,j,k)
						end do
					end do
				end do 
			end if 
		
			new_tracer%q(:,:,:,itracer) = rk_q
		  call fill_merid_halo_cell(new_tracer%q(:,:,:,itracer), all_halo=.true.)
			call fill_zonal_halo_cell(new_tracer%q(:,:,:,itracer), all_halo=.true.)
	  end do
		end associate

	end subroutine step_adv_ieva

end module step_mod